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We examine the crescent singularity of a developable cone in a setting similar to that studied by 
Cerda et al [Nature 401, 46 (1999)]. Stretching is localized in a core region near the pushing tip 
■ and bending dominates the outer region. Two types of stresses in the outer region are identified and 

shown to scale differently with the distance to the tip. Energies of the d-cone are estimated and the 
, conditions for the scaling of core region size R c are discussed. Tests of the pushing force equation and 

direct geometrical measurements provide numerical evidence that core size scales as R c ~ h 1,l3 R 2 ^ 3 , 
where h is the thickness of sheet and R is the supporting container radius, in agreement with the 
proposition of Cerda et al. We give arguments that this observed scaling law should not represent 
the asymptotic behavior. Other properties are also studied and tested numerically, consistent with 
' our analysis. 
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As we crumple a piece of paper in our hands, two types of singular structures appear in the crumpled paper: folding 
ridges and point-like vertices. Energies are condensed into a network of such singularities. The properties of ridges 
have been studied thoroughly. Scaling laws governing the energy and size of the ridge have been obtained analytically 
and tested numerically [6-9]. Point-like singularities are also studied extensively [1-4, 11, 12, 16, 17], however, current 
understanding of their properties is not as complete as that of ridges. 

In this paper, we consider a single conical vertex studied by Cerda et al [1, 2]. One experimental realization is to 
push the center of a circular elastic sheet of radius R p axially into a cylindrical container of radius R, as illustrated 
^1 in FIG. n This is the simplest volume-restricting deformation of the sheet and causes the center of the sheet to move 
into the container by a distance d. It is useful to express the deflection of the sheet by e = d/R. Due to the constraint 
of unstretchability, the sheet deforms into a non-axisymmetric conical surface which is only in partial contact with the 
^-j- ' edge of the container. In the limit that thickness h of the sheet goes to zero, since bending modulus (~ h 3 ) vanishes 
, faster than the stretching modulus (~ h), there would be pure bending over the sheet and Gaussian curvature would 
be zero everywhere. Mathematically, such a conical surface is called perfectly developable cone 01 (d-cone). In this 
limit, some models about the shape of the d-cone have been proposed [2-4, 17). These models only give outer region 
solutions of d-cone shape, in the sense that they do not consider the stretching energy that is inevitable on a real 
sheet with finite thickness. For a real sheet, it must stretch near the tip, because otherwise, the curvature at the tip 
would be divergent, since curvature goes as 1/r, where r is the distance to the tip, thus causing divergent energy. 
Therefore, it is the finite thickness that causes the sheet to stretch greatly in a small region near the tip. This small 
. region is called the core region. It is where energetically expensive stretching is localized and its size is governed by 
the competition of the bending and stretching energies of the whole cone. As a result of the stress focusing in the 
core region, crescent-like shapes come out where bending stresses are big, as shown in FIG. In addition, as we will 
O ■ see later, finite thickness also causes a small amount of strains in the outer region. 

One problem of great interest is concerned with the size of the core region, characterized by the radius of curvature of 
the crescents. We want to know whether there is scaling behavior of core size and to determine any scaling exponents. 
Intuitively, we expect that this size R c should have a dependence on thickness h, since R c goes to zero as h goes to 
zero. Cerda et al [lj propose that R c scales as R c ~ h^^R 2 ^. This says that besides h, the supporting container 
radius R also determines R c . This is surprising because stretching energies are supposed to be localized in the core 
region, so the core should not be able to know about the length of outside container radius. In this work we explore 
this problem as well as other properties of a d-cone. In Sec. II, we describe the energies and forces that give rise to the 
crescent singularity. We first study the elastic properties of a truncated d-cone formed by cutting a cap region at the 
center of a regular d-cone, and then investigate the energetic variations as we join the cap region into the truncated 
d-cone. From them, we discuss the conditions for the existence of scaling behavior of R c . Details of numerical models 
of simulating an elastic sheet and producing desired shapes are presented in Sec. III. After numerical results that 
agree with the predicted d-cone properties are shown, we use two different methods to look for scaling exponents of 
R c . Finally, the limitations of and implications from our findings and future work are discussed in Sec. IV. 



FIG. 1: A d-cone appears when we push the center of an originally flat sheet against the edge of a cylindrical container from 
below with force F. This is a typical simulated d-cone shape, with side length / = 60a, container radius R = 38a, displacement 
d = 0.157? and thickness h = 0.102a, where a is the lattice spacing. Two sets of coordinate systems are shown. 




FIG. 2: Sketch of the core region of a d-cone. Crescent area is shown in white. Directors are shown as black lines. Many 
directors converge toward the two tips of the crescent. Dashed circle tangent to the crescent at its center defines the radius R c . 



II. THEORIES 



We begin by stating the connection between the deformation of the sheet and its elastic energy. In a thin d-cone 
the dominant energy is due to the bending distortion in outer region far from the crescent. We discuss the form of 
this energy and the stresses associated with it. Next we outline the energetic effects of the crescent region. Finally, 
we focus on how each of these energies is influenced by a change in the crescent radius R c . We sketch how scaling 
properties of the energy produce scaling behavior in the sheet. 

We first specify two coordinate systems. We define x — y as the horizontal plane of the supporting cylindrical 
container's edge, with origin on the axis of cylinder and z axis pointing upward (see FIG.^I. We then define an x' — if 
plane as parallel to x — y plane but with origin at the tip of the sheet and z' axis pointing opposite to z axis, x' — y 
plane is displaced from the x — y plane by a distance d. Next we define {p, 9) as the polar coordinates in the x' — y' 
plane. In the cylindrical coordinates {p9z'}, d-cone is described by z' = pip (9), in the limit that thickness h — > 0. We 
shall denote this limiting surface as the asymptotic d-cone. 

In equilibrium, the sheet assumes a conformation that minimizes the elastic energy. Thus, the actual core radius is 
that which minimizes this energy. Two forms of energy must be considered, bending energy B and stretching energy 
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S. The bending energy density is proportional to the square of the total curvature C(r) 0]. This C(r) defined as the 
trace of the curvature tensor, is sometimes called the mean curvature |10| . The constant of proportionality is called 
the bending modulus n. Thus 

B=±nJdA C(rf , (1) 

where J dA denotes the integral over the surface. In general there is a second form of bending energy, proportional 
to the average Gaussian curvature. We show later in this section that this Gaussian curvature energy is unimportant 
for the present system. 

The stretching energy density is proportional to the square of the strain tensor 7. For an isotropic material, there 
are two forms of stretching energy with constants of proportionality Yh and Y\h ja| 

S=\J dA [Yh (Tr 7 ) 2 + Y x h (Det 7)] , (2) 

where Y is Young's modulus, and Y\ is a convenient combination of Y and dimensionless Poisson's ratio v. They are 
related to bending modulus through k = Yh 3 / 12(1 — v 2 ). 

The actual strain and curvature fields are those which minimize the B + S. The variational minimization amounts 
to a statement that the normal forces on each element must balance. This statement is known as the "force von 
Karman equation" 

d a dpM aP = a aj3 C a p + P , (3) 

where M a p are torques per unit length, a a p are in-plane stresses, and P is external force per unit area of the sheet. 

When we push a circular sheet of radius R p into a container to form a d-cone, strain and curvature distortions are 
created and elastic energy is stored. This energy arises from several effects. It is convenient to discuss it by creating 
the final surface in two stages. In the first stage, we assume a value for R c and cut a circular hole of that radius in 
the center of the sheet. We then force this perforated sheet into the container by exerting tension on the inner edge of 
the hole. In the unstretchable limit, h R, this shape becomes a truncated <i-cone, which serves as the outer region 
of a regular <i-cone. In this region, the curvature at every point vanishes in the radial direction. The total curvature 
C has the form C = 4>(9)/r. Its bending energy Bq is thus given by 
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Bq = t: k / rdr/r 2 I d6(/) 2 {9) (x Klog(R/R c ) . (4) 



In the limit /i < JJ C « i? p , this Bq dominates the energy, as we justify below. Then the function (f> takes the form that 
minimizes Bq subject to the constraints on the d-cone. In this geometry, the force von Karman equation simplifies 
greatly, and the form of 4>(9) as well as the associated stresses can be found explicitly 

We first determine the transverse stress agg using Eq. @. The nonzero components of torque tensor M follow from 
the constitutive law implicit in the energy equations [18j: M rr = nvCge — nv<f>{ff)/r and Mgg = nCgg = K(f>(9)/r. The 
force von Karman equation then reduces to d 2 M rr + d 2 Mgg/r 2 — aggCgg + P, from which we obtain 

K ( 6(6)\ Pr 

where the dots above functions denote 9 derivatives. In the same region where only Cgg is nonzero, there is no pressure 
from external force acting on the sheet, i.e. P — 0. Therefore, according to Eq. (JSJ, stress age goes as n/r 2 . These 
stresses arise from the requirement that the normal force due to changing torques is balanced by the normal force due 
to in-plane tension of the sheet. We call them type I stresses and denote them by symbol a^\ 

On the other hand, we consider the force balance of a region that encloses area between inner radius R c and outer 
radius r, where r can take values between R c and R. The tension exerted on the inner edge of this region is equivalent 
to the central pushing force of a regular <i-cone. This force must be balanced by the force due to radial in-plane stress 
a rr on the outer perimeter of the region. Let f3 be the angle between a radial line or generator of the d-cone and the 
horizontal. Since tan/3 = ip{9), we have sin/3 = ip(9)/ y^l + ip 2 (9). The balance of vertical forces yields 

a rr r sin (3d9 = F , (6) 



which holds for every value of r from R c to R. It is easy to see that type I stresses alone can't satisfy this equation, 
since type I stresses go as 1/r 2 , they would give a 1/r prefactor on the left side of equation, while the right-side of 
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equation is independent of r. Therefore, the integral from type I stresses must vanish and there must exist some 
additional stresses in the outer region that scale as 1/r to satisfy Eq. ©. We call these stresses as type II stresses 
and denote them by . They persist up to the supporting container edge, where normal force from the container 
counteracts the external pushing force. We can write cr^r — Fe(8), where e(9) is a function only of 9 and satisfies 
J e(9) sin (3(16 = 1. It is obvious that e(9) is of order unity. 

To estimate the magnitude of type II stresses, we notice that F = dE/dd = (dE/de)/ R « k/R. Thus u' 2 ' ~ 
F/r n/(rR). The type II stresses are comparable with type I stresses only near the container edge since the ratio 
(T ( 1 )y' cr ( 2 ) ~ r jJi f or R c < r < R. It is noted that they are both due to nonzero thicknesses. Since in-plane stress 
tensors are related to strain tensors through |fj a a/ 3 = [Yh/(l — v 2 )](j a i3 + ve ap e(3 T jp T ), we have strains 7«cr/ (Yh). 
The stretching energy of truncated (i-cone is then 

S ~ nh- 2 J ((<r (1) W 2) )/(17i)) 2 rdrdO w nh 2 / R 2 C . (7) 

We now examine the profile of external force along the container edge. Here the external normal force pressure 
P of Eq. © is nonvanishing. This normal force causes a small transverse deflection of the sheet and hence induces 
both curvature and strain near the edge. Noticing that the curvature induced is in radial direction and has opposite 
sign with Cgg, we expect mean curvature to be reduced near the edge. Due to the translational symmetry along the 
region of the surface in contact with the cylindrical container, the normal force pressure is found to be independent 
of azimuthal angle 9. However, besides this ^-independent term, Cerda and Mahadevan |0| show that a (5-function 
term of 9 emerges in the normal force pressure expression in order to make the torque balance, and this singularity 
happens at the take-off angular positions, where the sheet begins to bend away from the container. Since stresses and 
curvatures should avoid singular behavior for the sake of energy, according to Eq. (jSJ, there must exist a (5-function 
term in <j){9) to cancel that from normal force pressure P. Integrating over 9, we conclude that (f)(9) should have 
a jump at the take-off positions. This result is consistent with the geometrical requirement. More quantitatively, 
following Ref. ! we nn d that the ratio of the normal force contributed by the (5-function term to that contributed 
by the ^-independent term is tan# c /[2(7r — 9 C )}, where 9 C is the half aperture angle of non-contact region. For small 
deflections, taking 9 C 1.21 rad Q, we obtain the value of this ratio 0.69. We will numerically verify this ratio later. 

Having discussed stresses and forces, we now consider in detail the bending energy Bq. Bending energy comes 
from both total curvature and Gaussian curvature. Since the truncated ci-cone has no Gaussian curvature, we only 
need to take account of the total curvature contribution. The reduced curvature (f)(9) is related to the reduced height 
if)(9) = z'/p defined in the beginning of this section: (f)(9) = tp(6) + if>(9) - ip(9)ip 2 (9)/(l + i)j 2 (9)). For small to 
moderate deformations of the sheet (e < 0.4), C « (ip + if))/r ~ e/r. Throughout the paper we shall focus on this 
moderate range of e. The bending energy Bq has the form 

B ° = \ J ° 2<iA ~ ° lKe2 ln ,Rc) ' (8) 

where G\ is a geometrical factor. Comparing this with So from Eq. (JJJ), we see that bending dominates the outer 
region. 

Although Gaussian curvature is zero everywhere in a truncated d-cone, it is certainly not the case for a regular 
d-cone. The Gaussian curvature contribution to the bending energy is Bq = (kg/2) J KdA, where kq is Gaussian 
curvature coefficient, K is Gaussian curvature and A is the area. According to the Gauss-Bonnet theorem, the integral 
of Gaussian curvature over a region M of a surface is related to the integral of the geodesic curvature n g over the 
boundary of that region through J M KdA = 2tt — J dM n g ds. Choosing M to be the whole sheet without perforation, 
we get Bq — (kg/2)(2tt — J n g ds). If the sheet were a perfectly developable conical surface, the geodesic curvature at 
the perimeter would be the same as that of a regular cone of same size: n g = l/i? p , where R p is the perimeter-to-tip 
distance. However, due to stretching of the real sheet and the pushing of container edge, the shape is not perfectly 
developable so there are variations in the perimeter-to-tip distances along the boundary, which make n g deviate from 
1/Rp. Since stretching is small compared with bending and the sheet is nearly developable over most of the surface, 
this deviation is small compared with 1/R P , and it should depend on thickness h. We suppose that K g is a regular 
function of h, and may be expressed n g « (l/R p ) (1 + (hj R p )Q(9, e)), where 9 is a dimensionless function of 9 and e. 
Substituting this into the equation for bending energy, we obtain Bq = — (kq/2)(}i/ R p ) J 0(9, e)d9. Compared with 
Eq. ©, this energy is only a small fraction (h/R p )/ \n(R p /R c ) of the bending energy contributed by total curvature. 

The truncated d-cone discussed above has no crescent singularity. This singularity arises from the further constraints 
of filling in the hole in the truncated cone. The flat disk of radius R c that was removed must be distorted in order to 
join onto the truncated cone. This distorted cap exerts forces on the truncated (i-cone and distorts it in turn. All of 
these distortions add elastic energy to the energy Bq + So of the initial truncated d-cone. One obvious addition is the 
bending energy Bq due to Gaussian curvature discussed above. Another simple consequence of adding the core region 
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FIG. 3: (a) Sketch of a cross-section of the d-cone. Solid hypotenuse is for the d-cone when h — > and has slope e; dashed 
hypotenuse is for a real sheet (h 7^ 0) and has slope e' . (b) Sketch of fitting a hyperbola to find R c from geometry. 



is to alter the slope of the outer region generator from e to a higher value e' as shown in FIG. Eta). To relate e' to e, 
we imagine a cross-section of the d-cone extending from the container edge to the tip and back down to the edge (see 
FIG. Eta)). In the truncated <i-cone, the generators have the same slope e as the straight line from the container edge 
to the vertex, as shown by the solid lines in FIG. Eta). Now for the real sheet, the cap is curved, with a curvature 
radius of order R c , but still, its uppermost point remains at height d. This rounding at the top necessarily increases 
the slope of the outer cone generator, as illustrated by the dashed lines in the same figure. When deformation is 
small, e' s» d/(R — R c /2) ss e(l + R C /(2R)). The altered bending energy is given by Eq. JHJ with e replaced by e': 
B' « Gi«e 2 (l + R c / R) \n(R p /R c ). This increases the Bq energy by a factor (1 + R c /R). This addition to the energy 
favors small R c . 

Besides these variations, there are other forms of added energy. We denote them as AE. It consists of stretching 
energy AE added to the original region plus the energy E c of the cap or core region. The added energy AE is 
expected to be much smaller than the dominant energy Bo discussed above. Nevertheless, AE is all-important in 
determining the core-radius R c . 

Further distortions increase the energy of the outer region. The crescent singularity increases the curvature near 
the tips of the crescent. As shown in FIG. El the generators, initially distributed smoothly around the inner hole, are 
now concentrated at the tips of the crescent, increasing the curvature there. The generators splay outward from these 
tips. Observation suggests that a finite fraction of the generators are pushed to the tips. These generators appear to 
splay outward less rapidly when the size R is increased. These observations suggest that a) the crescent increases the 
exterior curvature energy, and b) this energy decreases as the size R — > 00 for a given R c . This energy penalty favors 
small R c /R. We do not have a more specific estimate for this energy at present. 

We now consider the energy E c of the cap or core region itself. In order to bridge the truncated cone, the 
average curvature C must be of order e/i? c . If the curvatures are uniform, the associated bending energy is of order 
k J dAC 2 ~ ne 2 . However, such uniform curvatures are not optimal. The smoothly-curved cap region has typical 
Gaussian curvatures of order e 2 /R 2 , consisting of a negative part in the buckled region and a comparable positive part 
elsewhere. The presence of Gaussian curvature induces strain of order e 2 , and a stretching energy of order Kh~ 2 e 4 R 2 . 
The system may alleviate this large stress energy by concentrating the curvature, as in the simple stretching ridge 
The observed crescent singularity is presumably the result of this concentrated curvature. The length of this crescent 
is of order R c ; its width w is evidently that which minimizes E c . If w is a fixed fraction of R c , the above estimates 
show that E c grows as R 2 . To improve this energy, w/R c must go to zero as R c grows. The total curvature is then 
of order 1/w, and the bending energy of the core goes as kwR c /w 2 ~ R c /w. This energy must grow indefinitely as 
R c — > 00. Evidently, the core energy favors small R c . 

With these energies in mind, we now investigate how R c is determined. It is clear that R c is controlled by the 
competition of energies. The total energy is E = Sq + B' Q + Bq + AE. Our purpose is to find an optimal R c that 
minimizes the total energy. First, taking derivative of B' Q with respect to R c , we have dB' /dR c — G\ne 2 [— 1/R C + 
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(ln(R. p /R c ) — 1)/R]. In the regime of interest, R is comparable with R p , and R c <§C R, so the first term dominates the 
second term: dB' /dR c w —GiK€ 2 /R c . This implies that £? favors large i? c . Next, we notice that since OBq/DRc 
and dSo/dR c are much smaller than dB' Q /dR c in this regime, the competition happens mainly between outer bending 
energy B' Q and the added energy AE. Though we know little about the form of the added energy AE, the above 
observations suggest that a) it favors small R c , b) it decreases as R increases, and c) it involves stretching, and thus 
increases with decreasing thickness h. To illustrate the possible effect of this energy, we posit that AE is homogeneous 
in R c , R, and h 

AE = Kg(e)R s c iT* h~ s+t , (9) 

where g{e) is an unknown function of e and s and t are unknown positive exponents. The h power law is determined 
from s and t by dimensional consistency. The energy minimization with respect to the variational parameter R c then 
implies 

BF 

= i? C 7r— = —G\ne 2 + sAE . (10) 

Oti c 

Using the assumed power-law form of AE, we infer 

R c ~ Rt/sh 1 -*/ 3 . (11) 

Thus the assumption of Eq. © implies that R c should increase as a power of R. Since R c cannot exceed R, the 
power must be smaller than unity, so that t must be smaller than s. 

Though this assumed form for AE gives a pleasing result, we warn that it cannot be correct in the limit of large 
R/h. If it were correct, it would imply that R c /h goes to infinity with R/h. This implies that the core energy E c 
must go to infinity relative to ne 2 , as shown above. The increase of E c occurs irrespective of R. Thus E c must grow 
to dominate the total energy. Such a large E c is not compatible with a minimum of total energy. It must be balanced 
by some other energy favoring large R c . No such energy is apparent. Thus R c cannot grow indefinitely. Conversely, 
any observed power-law growth of R c must cease for sufficiently large R/h. In the numerical section, we will give an 
empirical discussion of the range of R/h in which the scaling of R c holds. 

These conclusions contrast with the scaling argument proposed in the initial studies of the d-cone Q 117| . Thus it is 
important to revisit this scaling argument and show how the contradiction with our result arises. The authors attribute 
the scaling to an energy balance between bending and stretching in the core region. To estimate the stretching energy, 
they consider a director traversing the sheet from one edge, through the core and on to the opposite edge, such as the 
dashed line in FIG.|3Ja). They note that increased R c increases the length of this chord (relative to the solid line) 
by a fraction of order (R c /R) 2 . If this strain is presumed uniform, the associated stretching energy in the core is of 
order nh~ 2 R 2 (R C / R) 4 . They balance this energy against the bending energy in the core, presumed to be of order 
ne 2 . This balance yields R c ~ R 2 / 3 ^/ 3 . 

This argument seems questionable, since it ignores energies much larger than the ones it includes. Its assumption 
of uniform strain leads to an stretching energy of uniform density outside as well as within the core. The argument 
includes that part of this energy lying within the core while ignoring the much larger part outside the core. If one tries 
to repair the estimate by supposing the excess length resides only in the core, then the strain becomes independent 
of R and the R c must be independent of R. In addition the argument ignores the large stretching energy in the core 
arising from the necessary Gaussian curvature within the core itself, as discussed above. Likewise, it takes no account 
of the strong difference between the two principal curvatures in the crescent. 

Our analysis of the energy given above allows us to infer the scaling of R c based on measurements of the central 
pushing force, F. The inference is valid so long as the dominant energy is the Bq energy above. The energy balance 
equation l|10|) allows us to write AE = G^kc 2 , where G2 is a numerical constant. Thus 

E « G\ne 2 ln{R p /R c ) + G 2 ne 2 . (12) 

The pushing force follows from F = dE/dd = (dE/de)/R. To explore the scaling of R c , we write R c as R c = 
h p R q RpA(e), where p, q and A are scaling exponents to be determined, and A(e) is a function only of e. Dimensional 
consistency requires that p + q + A = 1 . Previous work Q did not consider the possibility that R c has a dependence 
on R p . Here we address it by introducing A. From Eq. I|12l) and the scaling expression of R c , the pushing force is 
calculated to be 

F^^^(-p]nh-qlnR+(l-X)lnR p + f(e)) , (13) 

where /(e) = — \nA(e) + G2/G1 — eA(e)/2A(e) is a function only of e. We will use this equation to test the scaling 
relations in the numerical section, which follows below. 
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III. NUMERICS AND FINDINGS 

A. Numerical Model 

We model an elastic sheet by a triangular lattice of springs of un-stretched length a and spring constant k after 
Seung and Nelson 10] . Bending rigidity is introduced by assigning an energy of J(l — ri\ ■ rix) to every pair of 
adjacent triangles with normals n\ and n-x- When strains are small compared to unity and radii of curvature are large 
compared to the lattice spacing a, this model is equivalent to an elastic sheet of thickness h = aycBJ/fc made of an 
isotropic, homogeneous material with bending modulus k — J^/3/2, Young's modulus Y = 2ka/h^/3 and Poisson's 
ratio v = 1/3. Lattice spacing a is set to be 1. The shape of the sheet in our simulation is a regular hexagon of side 
length R p . The typical value of R p is 60a. 

To obtain a single d-cone shape, we need to simulate the constraining container edge and pushing force. As shown 
in FIG.m the edge lies in the x — y plane and is described by equation x 2 + y 2 = R 2 . Pushing in the center of the sheet 
is accomplished by introducing a repulsive potential of the form U{ OIce (zi) — —Fzi, where z\ is the z coordinate of the 
lattice point in the center and F is the magnitude of the pushing force. This force is applied in the p ositive z direction. 
The constraining edge is implemented by a potential of the form C/ e d g c = C p H(zi) / (((^J x 2 + y 2 — R) 2 + z 2 ) 4 + £ 8 ), 
where £, C p are constants and the summation is over all lattice pionts with coordinates (xi,yi, z{). H(z) is the unit 
step function, which makes certain that this potential only acts on the lattice points that have already moved into the 
container (those with z t > 0). The force associated with it decays rapidly once the lattice points go away from the 
edge. The conjugate gradient algorithm jl4j is used to minimize the total elastic and potential energy of the system 
as a function of the coordinates of all lattice points. FIG.n snows one such minimized configuration of the lattice 
grid. The simulated sheet is indeed only in partial contact with the container edge. In the rest of this section, we 
shall first compare our results to the d-cone predictions of previous work. The good agreement indicates that our 
numerical realization is reliable. Next we shall investigate the scaling behaviour of R c . 

B. Curvature Profile 

We first test our simulation by measuring the curvature profile of the sheet, and comparing it to the prediction 
under certain limiting situations. We determine the curvatures approximately from each triangle in the sheet. For 
this measurement, we take the curvature tensor to be constant across each triangle. We calculate it using the relative 
heights of the six vertices of the three triangles that share sides with the given triangle 0] . The six relative heights 
Wi normal to the triangle surface are fit to a function of the form 

Wi = h + b 2 u t + b 3 Vi + b±u 2 + b^Vi + b 6 v 2 , 2 = 1, — ,6 (14) 

where {v,i,Vi, u>i} are coordinates of the vertices in a local coordinate system that has w axis perpendicular to the 
surface of the given triangle. This choice of local coordinate system ensures that 62 and 03 are negligible so that 
curvature tensors can be determined only from the coefficients of quadratic terms. In practice, our numerical findings 
do show that the values of 02 and 63 are on the order of 10~ 2 or lower. Therefore, curvature tensors follow immediately 
from the identification C uu = 2 x 64, C vv = 2 x & 6 , C uv = 65. Total curvature C is defined as the trace of curvature 
tensor: C = C uu + C vv . FIG. Ola) gives a typical plot of total curvature versus azimuthal angle on a d-cone surface at 
three different distances from the tip. Curvature is measured in units of 1/r and taken to be negative in the convex 
region, where the sheet touches the supporting container. The three curves mostly collapse into one single curve. This 
roughly verifies that curvature goes as 1/r. However, near 9 = 0, the normalized curvature becomes systematically 
smaller for smaller r. We attribute this departure to the influence of the stretched crescent region. 

The shapes of the curves are in qualitative agreement with our analysis below. The surface is in contact with the 
container for angles 9 larger than some 9 C in magnitude. For \9\ > 9 C , the normalized curvature at all r shown is 
constant and its value is that of a simple cone of the same e. As \9\ becomes less than the "take-off" angle 9 C , the 
surface curves inwardly away from the container and the total curvature begins to increase in negative direction. The 
crescent singularity is most likely to live near the region where total curvature reaches its negative maximum. As \9\ 
decreases further, the curvature goes to zero, changes signs, and reaches a central positive maximum at 9 = 0, which 
corresponds to where the sheet is maximally deflected. 

As thickness h and deformation e go to zero, Cerda and Mahadevan obtained the exact solution of curvature profile 
0^3- Although it is impossible for us to get the h—*0 profile, we want to compare our data with their prediction 
in the limit that e — 0. To do so, recall that total curvature, as mentioned in Sec. I, is related to the shape through 
C = [Tp{9)+'4){9)-ij(9)^ 2 (9) / {\ + ^ 2 {9))]/r. Since ip(6) cx e, we have Cr/e = a Q + a 1 e 2 + a 2 e 4 + 0(e 6 ), where a ,ai and 
02 could depend on 9. Hence, from this equation, we can extrapolate the curvature at e = by fitting the curvature at 
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three known values of e with basis functions {1, e 2 , e 4 }. In FIG.^fb) we plot curvature profiles at e = 0.20, 0.15, 0.10 
and 0. The e — *■ profile is extrapolated from other three profiles using the fitting theme discussed above. Solid line 
is the theory curve of exact solution. As we can see from this graph, as e decreases, the curvature profile becomes 
closer to the theory curve that is based on the assumption that e = and h = 0. Our e = profile agrees well with 
the theory curve in the central peak region. However, striking difference still exists in the range between the take-off 
positions and negative maxima. This discrepancy can be explained by the effect of nonzero thickness. For real sheet, 
h =/= 0, from Eq. (J5J, we observe that <f> tends to avoid jump at take-off positions in the non-contact region where 
external pressure P = 0, otherwise there would be singularity in the strains. Therefore the observed curvature profile 
has to be rounded rather than having an abrupt kink in slope, which requires that curvature profile be curved up as 
the sheet leaves the container, just like what is displayed in the plot. Besides, the same reason can explain the slight 
difference near take-off points in the curvature profiles at three distances in FIG.^Ja). For smaller r, the curvature 
may be more influenced by nonzero thickness effect, so the profile at r = 30a is more rounded than the profile at 
r = 50a. In addition, we note that the theory curve based on the model proposed by Chaieb et al 4| does not match 
our data. 
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FIG. 4: (a) Azimuthal profile of normalized total curvature on a d-cone at three fixed distances from the tip. Curvature is 
positive for the concave region; negative for the convex region. Curvature is normalized by 1/r, where r is the fixed distance 
from the tip. This plot is for the d-cone shape shown in FIG.^at distances 30a (circle), 40a (plus) and 50a (square) from the 
tip. The container touches the d-cone surface at r = Ry/1 + e 2 £s 38.4a. (b) Normalized curvature profiles for four different 
values of e at fixed thickness h = 0.102a. Curvature is normalized by e/r. Profile at e — > is obtained by extrapolating from 
profiles at e = 0.10, 0.15 and 0.20. Solid line is the exact solution as both h — *• and e — ► [T^l . 

From the curvature file, we can measure the opening angle of d-cone. The opening angle is defined as the angular 
distance between the two "take-off points" where the sheet loses contact with the container. For asymptotically thin 
sheets with small deformation, this angle is predicted to be 138° from the exact solution, as illustrated by the solid 
line of FIG. Bib). As e decreases to zero, the opening angle measured from curvature profiles tends to converge to 
that indicated by the solid line. We expect the agreement with the prediction when the elastic thickness h vanishes. 
This value was approximately confirmed by experiments, which yielded 130° |2j. 

C. Normal Force 

The azimuthal profile of the normal force pressure from the container is displayed in FIG. |S] It is evident that 
we do observe sharp peaks at the take-off positions, which confirm the r5-function term in the normal force pressure 
proposed in 01 • The pressure drops quickly to zero as one enters the buckled region where the sheet bends away from 
container. The angular separation of the take-off positions is about 2.21 rad = 127°. To verify that the sharp peaks 
in the normal force pressure have the proper strength, we calculate the ratio of the normal force from the (5-function 
term to that from the ^-independent term. The ratio is found to be 0.70 from our data, compared well with the 
theoretical prediction 0.69 obtained in Sec. II. 
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Our measurements in Sec. III.B and this section confirm that our simulation accurately represents an elastic sheet 
as desired. We now report the observed behavior of the core region. 




-3-2-1 1 2 3 

e (rad) 



FIG. 5: Azimuthal profile of normal force pressure. This is for the shape shown in FIG. Normal force pressure is measured 
in relative units. The right peak is located at 1.12 rad; the left peak is located at —1.09 rad. 



D. Scaling of Core Size 

To study the scaling law of core region size, it is natural to start estimating core size by finding the radius of 
curvature of crescents. From azimuthal profiles of curvature at different distances, we find the locations of the triangles 
with maximal negative curvature on both sides at each fixed distance. Centers of these triangles are projected onto 
x — y plane. The best fitted straight lines to the projection points on two sides serve as the asymptotic lines of a 
hyperbola, and the central forcing point is taken as the midpoint of the same hyperbola. This process is illustrated in 
FIG.EJb). R c is the radius of curvature at the midpoint of the fitted hyperbola shape. FIG.Elshows the dependence 
of R c measured in this way on h and R when deformation e is fixed at 0.10 and 0.15. We can observe power law 
dependence from these plots. The power law fits give 0.334 ± 0.022 and 0.380 ± 0.014 for thickness dependence, and 
give 0.574±0.021 and 0.597±0.032 for radius dependence. These values roughly agree with the 1/3 and 2/3 exponents 
proposed in Ref. 0. 



(a) (b) 




FIG. 6: The plots of (a) R c versus h and (b) R c versus R at e = 0.10 and e = 0.15. R is fixed for plots in (a); h is fixed for plots 
in (b). The straight lines are power law fits. The fitted values of power in (a) are 0.334 ± 0.022 for e = 0.10 and 0.380 ± 0.014 
for e = 0.15. The fitted values of power in (b) are 0.574 ± 0.021 for e = 0.10 and 0.597 ± 0.032 for e = 0.15. The error bars in 
the graph come from the numerical hysteresis effect when thickness of sheet or radius of container is increased and decreased 
through the same value during the energy minimization process. 
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As pointed out above, however, the lattice model can only accurately simulate an elastic sheet where the radius 
of curvature, 1/C, is locally much greater than the lattice spacing. We expect this condition not to be well satisfied 
in the core region, where singularity happens. Indeed, our data indicate that Ca can be as big as 0.5 in this region 
even for small e. Thus we have no hope of maintaining complete accuracy in this region. In addition, we find that 
the curvature data do not exhibit the kind of shape as shown in FIG. 0] at small distances. Therefore, in the above 
procedure of measuring R c geometrically, we have to use curvature profiles at big distances, which makes it not so 
accurate to determine R c in this way since R c is supposed to be determined from the information in the neighbourhood 
of the tip. 

Our second approach of looking for the scaling relations is to use the force equation l|13|) . In our simulation, we 
keep the spring constant k fixed (k = 1), which is equivalent to fixing two-dimensional Young's modulus (Yh) since 
we have Yh = 2fca/v / 3- Hence k cx Yh 3 cx h 2 . From Eq. I|13fl we have 



h 2 

F ^ —(-ph-ih-q\nR+{p + q)h-iR p + D) , (15) 
rslx 



where B and D only depend on e, which we now fix at 0.10. Fixing e is realized through a process of several 
minimizations. We slowly adjust the pushing force until the desired e value is reached. In every step of the process, 
the previously obtained minimized configuration is used as the input for the next minimization procedure. 

To find the relations between exponents, we first fix h by fixing the bending coefficient J defined in Sec. III. A, and 
measure force F on the minimized shapes for different values of R. FIG. 0a) gives the semilog plot of FR versus R 
when h is fixed at 0.102a and R p is fixed at 60a. The linear feature of the plot agrees with the prediction from Eq. 
(|15|l . We denote the slope and intercept of the best fitted line by Si and I±, respectively. Then Si = — 0.0342±0.0011 
and h = 0.1714 ± 0.0042. From Eq. O we find 

hB = h 2 (-pin h +(p + q) In R p + D) (16a) 
SiB = h 2 {-q) (16b) 

Next, for the same value of e, R = 37a and R' p = 60a are fixed while h is varied as the corresponding equilibrium 
force F is measured. F/h 2 versus h semilog plot is shown in FIG.EJb). Similarly the slope of the best fitted line 
S 2 = -0.0494 ± 0.0005 and intercept I 2 = 0.0108 ± 0.0006, which according to Eq. (JSJ are given by 

hB = -(-qlnR+(p + q)lnR'+D) (17a) 
R 

S2B = (17b) 

Notice that Eqs. (|16afl . (|16bl) . 1)1 7al) and l|17b|) constitute a homogeneous system of four linear equations in four 
variables {p, q, D, B}. In order to have non-trivial solutions, the determinant of its coefficient matrix must vanish. 
To check this condition, we write out the determinant of dimcnsionlcss matrix of coefficients explicitly 



(In Rp - In h) 


\nR' p 
1 



In Rp 1 -h/ti 2 
1 S1//1 2 

(InR'-lnR) 1 -I 2 R 







S 2 R 



(lnR/h 2 )Sx + h/h 2 - [R\n{hR' p / R p ))S 2 - Rh 



which is calculated to be 0.016±0.861, indeed including 0. In addition, we obtain the relations between these variables: 
q = (1.879 ± 0.428)p and B = (0.547 ± 0.005)p. 

The third relationship we can extract from force equation l|15|l is F versus R p . To test this, we keep the hi = 
0.1549a and R\ = 25a fixed while measuring pushing forces for different R p . The plot is shown in FIG. Efc). 
On the one hand, the slope of the best fitted line is 0.00461 ± 0.00003 and intercept is -0.0147 ± 0.0001. On 
the other hand, we can calculate the slope and intercept of this graph using the relations between p, g, D and B 
obtained above. Specifically, from Eq. i|15|) we have slope S3 = h 2 (p + q)/{BRi) = 0.0051 ± 0.0008, and intercept 
^3 = h\/ {BRx){—phxhx — glni?i + D) = —0.0157 ± 0.0026, both in good agreement with the direct fitted values 
from plot. These tests verify the self-consistency of our data and the validity of the form of the force equation, thus 
supporting the scaling behaviour of R c . 

The useful information we obtain from these tests concerning the scaling exponents is the relation between p and 
q. We can not determine the value of A in Eq. i|13|) just from these tests. However, geometrical measurements show 
that R c varies very little with R p . Indeed, we find that the standard deviation of values of R c at different R p is only 
about a couple percents of the mean value of R c . It seems this observation serves as a good footing for us to believe 
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that R c is independent of R p . If we make this assumption, that is, if we take A = 0, then from the relation between 
p and q, we can obtain p — 0.355 ± 0.053. This is close to the corresponding value from geometrical measurements. 

Besides e = 0.10, we perform the above tests of the force equation for other values of e. The results are summarized 
in TABLE[I] We notice that all the values of p are consistent with scaling exponent of 1/3. In addition, by comparison 
of Eq. (|13fl and Eq. i|15|) . it is easy to see that B is inversely proportional to e. This relation can be readily confirmed 
from the data shown in the table. 

As mentioned in Sec. II, R c scaling is expected to break for sufficiently large R/h. However, in our simulations 
discussed above, all the data are consistent with the assumption that there exists a scaling law of R c , which indicates 
that the asymptotic regime of large R/h is not yet reached. The values of R/h in our simulations vary between 100 
and 500. We conclude that the lower limit of this asymptotic regime should be at least above 500. It is unsettling 
that the asymptotic regime should require such large R/h, but similar behavior occurs in relatedphenomena. For 
example, the asymptotic scaling of the stretching ridge requires a ratio of R/h of over a thousand |9(. 

The reason of this unusually high lower limit may be the following. The minor radius of the crescent is expected to 
be asymptotically much smaller than R c . However, our observation shows that this minor radius is comparable with 
R c for the range of R c / h we are able to simulate. In this respect we see directly that our system is not asymptotic. We 
have noted in Sec. II that the system's asymptotic behavior depends strongly on how this minor radius behaves. Thus 
one cannot expect asymptotic scaling of R c without asymptotic behavior of the minor radius. Since our simulations 
do not reach this behavior, we should not be surprised that R c may not show the asymptotic scaling, either. To 
demonstrate that the minor radius is not asymptotic, we show its behavior in FIG.|S] This plot shows the intersection 
of the surface with a vertical plane passing transverse to the crescent. Thus the radius of curvature of this curve at its 
peak is the minor radius. Both the horizontal and vertical scales are normalized by R c . It is evident from the plot that 
this normalized width is on the order of unity for all of the four thicknesses; this means the minor radius of crescent 
is comparable with R c . In addition, the feature that the minor radius of crescent relative to R c is growing smaller as 
thickness decreases provides evidence that the observed scaling behavior of R c in our simulations is non-asymptotic. 




FIG. 7: Force plots when e is fixed at 0.10. (a) FR versus R semilog plot. The data is fitted to be FR = -0.0342 x In i?±0.1714. 
The parameters we use are k = 1, J = 0.0013, R p = 60a. (b) F/h 2 versus h semilog plot. The data is fitted to be 
F/h 2 = -0.0494 x In h + 0.0108. The parameters we use are k = 1, R = 37a, R p = 60a. (c) F versus R p semilog plot. The 
data is fitted to be F — 0.0046 x ln_R p — 0.0147. The parameters we use are k = 1, J = 0.003, R = 25a. The relation between 
exponents p and g is g = 1.879p. 



TABLE I: Results of testing force equation for different values of e 



e 


Determinant 


B/p 


P 


0.10 


0.016 ±0.861 


0.547 ±0.005 


0.355 ± 0.053 


0.15 


-0.053 ± 0.973 


0.367 ±0.003 


0.344 ± 0.040 


0.20 


0.056 ± 2.076 


0.267 ±0.001 


0.381 ± 0.071 


0.25 


0.163 ±2.720 


0.203 ±0.001 


0.410 ± 0.073 
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FIG. 8: Intersections of a d-cone with a vertical plane consisting of its maximally deflected line and z-axis, for four different 
thicknesses. Only the region near the peak is shown here. Both scales are normalized by the corresponding R c for values of h 
as indicated in the legend. Thicknesses are in the units of lattice spacing a. 



IV. DISCUSSION 



In this paper we have explored properties of the conical singularity of a developable cone, especially the scaling of 
core region size. We have found out two types of stresses in the outer region of a single d-cone that scale differently 
with the distance to the tip. One type of stress arises from the normal force balance of the sheet without the external 
load, and scales as 1/r 2 . The other type of stress scaling as 1/r is needed to balance the external pushing force. 
However, it is revealed that the contribution of both stresses to the stretching energy is negligible compared with 
bending energy of the outer region. 

We also examined the normal force pressure from the container. The jump in the curvature derivative <fi/r requires 
a (^-function pressure from the container edge at the take-off points. We verified that this singular pressure is present 
with the predicted magnitude [Tt| |. As a consequence, a substantial fraction of the container forces comes from this 
(5-function term. 

Our numerical tests of pushing force equation suggest the existence of scaling behavior of R c in the regime we 
studied, that is, when R is comparable with R p . We obtain a simple proportionality factor between exponents p 
and q. Since R c has dependence on h, i.e. p is nonzero, it follows that q is nonzero, which implies R c must have 
a dependence on R. This is somewhat counter-intuitive. Moreover, geometrical measurements of core size provide 
suggestive numerical evidence that R c is independent of R p . By taking this assumption, we are led to a scaling law 
suggesting R c ~ ti^^R 2 ^ . Although our numerical results are consistent with the scaling prediction of Ref. Q, we 
were unable to justify the arguments leading to their prediction. The factors determining R c are necessarily subtle, 
since the dominant energy, Bq, depends only logarithmically on R c . Until a clear justification for this scaling can be 
found and validated, the apparent scaling we observed must be viewed as provisional. 

Our work in progress aims to explore the energetics leading to R c scaling in more details. Our preliminary findings 
suggest new features that may help to resolve this issue. First, one may construct variants of a d-cone that require no 
forcing at the core. Our preliminary data show that these variants have qualitatively different scaling behaviour than 
the conventional d-cones. Second, conventional d-cones appear to obey an unanticipated constraint at the container 
edge. The mean curvature appears to vanish there for a wide range of d-cone shapes. We anticipate that the energy 
focusing in d-cones will prove to be rich and revealing. 
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